###########################################################
### Analysis: Estimated Values Full Data #################
#########################################################

## Purpose: This script includes the models for the
## expected values for the full data analysis (Figures 2 and
## 3 in the main text). 
## Data In: 
## 1) data_files/datafull_5_6_2025.rds
## Data Out:
## 1) figures/race_time.pdf
## 2) figures/cohort_time.pdf
## 3) figures/party_time.pdf
## 4) figures/ses_time.pdf
## 5) figures/ses_party_intx
## 5) figures/ses_party_intx_bw
## Software Versions:

## R version 4.4.1

library(RColorBrewer) ## 1.1.3
library(tidyverse) ## 2.0.0
library(viridis) ## 0.6.5
library(MASS) ## 7.3.60.2
library(ggthemes) ## 5.1.0


## set working directory to 
## replication file

data.full <- readRDS("data_files/datafull_5_6_2025.rds")


'%ni%' <- Negate("%in%")

#########################################
### Covariate-Year Interaction Models ###
#########################################

data.cohortintx <- data.full %>% filter(cohort2 != "Before Greatest" &
                                          cohort2 != "Millenials" &
                                          year >= 1987)


full1 <- lm(inequality.num ~ year_5*ses, data = data.full)
full1a <- lm(inequality.num ~ year_5*ses + party.harm + cohort2 + race2 + gender, data = data.full)
full2 <- lm(inequality.num ~ year_5*party.harm, data = data.full)
full2a <- lm(inequality.num ~ year_5*party.harm  + ses + cohort2 + race2 + gender, data = data.full)
full3 <- lm(inequality.num ~ year_5*cohort2, data = data.cohortintx)
full3a <- lm(inequality.num ~ year_5*cohort2  + party.harm + ses + race2 + gender, data = data.cohortintx)
full4 <- lm(inequality.num ~ year_5*race2, data = data.full)
full4a <- lm(inequality.num ~ year_5*race2 + party.harm + ses + cohort2 + gender, data = data.full)


####################################
### Creating Small Multiples Plot ##
### of Estimated Values for each ###
### model ##########################
####################################

## We create Figure 2 in the main
## text in this block of code. Figure 2
## examines the percent of Americans who 
## agree the rich are getting richer by the 
## levels of our main predictor variables
## and time. 

## this function calculates s.e. and point estimates
## across simulation draws, used below

estimates.fct <- function(frame){
  est <- as.data.frame(matrix(NA, ncol = 4, nrow = ncol(frame)))
  colnames(est) <- c("Est", "Lower", "Upper", "Model")
  est$Est <- base::apply(frame, 2, function(x){mean(x)})
  est$Lower <- base::apply(frame, 2, function(x){quantile(x, .025)})
  est$Upper <- base::apply(frame, 2, function(x){quantile(x, .975)})
  est$Model <- colnames(frame)  
  return(est)
}

## in the code below and in all
## analyses in this script we calculate 
## estimated the values: average values of our
## outcome variable at different levels of our
## predictor variables, holding all other variables
## constant.
## We use our estimated regression coefficients to calculate
## these estimated values. 
## We estimate uncertainty around
## the estimated values using a parametric bootstrap. 

### model 1: year*ses
### simulating draws from the estimated
## distribution of model coefficients, centered
## at the observed estimates
vcov <- vcov(full1a)
betas <- full1a$coefficients

set.seed(08540)
beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)

## we create matrix to allow us to calculate estimated
## values for each set of simulated coefficients
## number of rows equal to
## number of estimates: levels of ses * levels of year
length(unique(data.full$year_5))
length(unique(data.full$ses))

matrix <- matrix(0, nrow = 50, 
                 ncol = length(beta.draws[1,]))
matrix[,1] <- 1
colnames(matrix) <- names(beta.draws[1,])
rownames(matrix) <- 1:nrow(matrix)

## year - different year for each level
## main year effects:

## for each level:
## each level gets a row of zeros (for first 1966-1970)
## and then a single year 
for(i in 2:10){
  matrix[i, i] <- 1
  matrix[i+10, i] <- 1
  matrix[i+20, i] <- 1
  matrix[i+30, i] <- 1
  matrix[i+40, i] <- 1
} 

## year*ses interaction effects:

for(i in 2:10){
  matrix[i+10, i + 22] <- 1
  matrix[i+20, i + 31] <- 1
  matrix[i+30, i + 40] <- 1
  matrix[i+40, i + 49] <- 1
}

## main ses effects:
matrix[11:20, "ses2"] <- 1
matrix[21:30, "ses3"] <- 1
matrix[31:40, "ses4"] <- 1
matrix[41:50, "ses5"] <- 1

## setting other variables at same leveL:
## Republican 
## White
## Boomers
## gender
matrix[, "race2Yes"] <- 1
matrix[, "gendermale"] <- 1


## for each draw of coefficients from the
## estimated distribution, we calculate
## estimated values using our matrix

sim <- base::apply(beta.draws, 1, function(x){
  vector <- matrix %*% x
  return(vector)
})
sim <- t(sim)
sim <- as.data.frame(sim)



## calculate 2.5th and 97.5th percentile
## of estimated values for confidence intervals
estimates.ses <- estimates.fct(sim)
estimates.ses$frame <- "ses"

estimates.ses$year <- rep(levels(data.full$year_5), 5)
estimates.ses$year_date <- gsub("-[0-9]{4}",
                                "", 
                                estimates.ses$year)
estimates.ses$year_date <- as.Date(paste0(estimates.ses$year_date,
                                          "-01-01",
                                          sep = ""))


estimates.ses$value <- c(rep(levels(data.full$ses)[1], 10), 
                         rep(levels(data.full$ses)[2], 10),
                         rep(levels(data.full$ses)[3], 10),
                         rep(levels(data.full$ses)[4], 10),
                         rep(levels(data.full$ses)[5], 10))

estimates.ses$year_date <- as.Date(estimates.ses$year_date)


ggplot(estimates.ses,
       aes(x = year_date,
           y = Est,
           color = value,
           group = value,
           ymin = Lower,
           ymax = Upper)) +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  geom_line(position = position_dodge(2 * 365)) +
  scale_color_viridis_d() +
  scale_y_continuous(limits = c(0, 1),
                     labels = scales::percent_format()) +
  theme_bw() +
  labs(x = NULL,
       y = "Estimated Percnet of Respondents\nWho Think 'Rich Getting Richer'",
       color = "SES Bin")
## saved as figures/ses_time.pdf
## 5x7
## included as top left panel of figure 2
## for caption: other variables held
## constant: Republican, Boomer, White, Male
## points slightly offset to aid in inspection

## estimates for results section 4.1
## 1991-1995
estimates.ses$Est[estimates.ses$year == "1991-1995" &
                    estimates.ses$value == 1] -
  estimates.ses$Est[estimates.ses$year == "1991-1995" &
                estimates.ses$value == 4]
## 12 pp difference

## non-parametric estimate for robustness:
data.full %>%
  filter(year_5 == "1991-1995") %>%
  filter(ses %in% c(1,4)) %>%
  group_by(ses) %>%
  summarize(mean = mean(inequality.num))
.881 - .742 ## 13.9 p.p. difference


## 1996 - 2000
estimates.ses$Est[estimates.ses$year == "1996-2000" &
                    estimates.ses$value == 1] -
  estimates.ses$Est[estimates.ses$year == "1996-2000" &
                      estimates.ses$value == 4]
## 13 pp difference

## non-parametric estimate for robustness:
data.full %>%
  filter(year_5 == "1996-2000") %>%
  filter(ses %in% c(1,4)) %>%
  group_by(ses) %>%
  summarize(mean = mean(inequality.num))
.844 - .681 ## 16.3 p.p difference


estimates.ses$Est[estimates.ses$year == "2001-2005" &
                    estimates.ses$value == 1] -
  estimates.ses$Est[estimates.ses$year == "2001-2005" &
                      estimates.ses$value == 4]
## 16.1 pp difference 

data.full %>%
  filter(year_5 == "2001-2005") %>%
  filter(ses %in% c(1,4)) %>%
  group_by(ses) %>%
  summarize(mean = mean(inequality.num))
.821 - .626 ## 19.5


## model 2: year*party
vcov <- vcov(full2a)
betas <- full2a$coefficients

set.seed(08540)
beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)

## create matrix for estimated values
matrix <- matrix(0, nrow = 30, 
                 ncol = length(beta.draws[1,]))
matrix[,1] <- 1
colnames(matrix) <- names(beta.draws[1,])
rownames(matrix) <- 1:nrow(matrix)

## year - different year for each level
## main year effects:

## for each level:
## each level gets a row of zeros (for first 1966-1970)
## and then a single year 

for(i in 2:10){
  matrix[i, i] <- 1
  matrix[i+10, i] <- 1
  matrix[i+20, i] <- 1
} 

## year*party interaction effects:
for(i in 2:10){
  matrix[i+10, i + 22] <- 1
  matrix[i+20, i + 31] <- 1
}


## main party effects:
matrix[11:20, "party.harmDemocrat"] <- 1
matrix[21:30, "party.harmIndep/Other"] <- 1

## other variables
## set at Boomers
## white
## ses: bin 3 (not ordered factor)
## gender - male
matrix[,"race2Yes"] <- 1
matrix[,"ses3"] <- 1
matrix[, "gendermale"] <- 1


## for each draw of coefficients from the
## estimated distribution, we calculate
## estimated values 

sim <- base::apply(beta.draws, 1, function(x){
  vector <- matrix %*% x
  return(vector)
})
sim <- t(sim)

## calculate 2.5th and 97.5th percentile
## of estimated values for confidence intervals
estimates.party <- estimates.fct(sim)
estimates.party$frame <- "party"

estimates.party$year <- rep(levels(data.full$year_5), 3)
estimates.party$value <- c(rep(levels(data.full$party.harm)[1], 10), 
                           rep(levels(data.full$party.harm)[2], 10),
                           rep(levels(data.full$party.harm)[3], 10))


estimates.party$year_date <- gsub("-[0-9]{4}",
                                  "", 
                                  estimates.party$year)
estimates.party$year_date <- as.Date(paste0(estimates.party$year_date,
                                            "-01-01",
                                            sep = ""))

estimates.party$year_date <- as.Date(estimates.party$year_date)

ggplot(estimates.party,
       aes(x = year_date,
           y = Est,
          shape = value,
           color = value,
           group = value,
           ymin = Lower,
           ymax = Upper)) +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  geom_line(position = position_dodge(2 * 365)) +
  scale_color_manual(values = c("#2c7fb8",
                                "#31a354",
                                "#f03b20")) +
  scale_shape_manual(values = c(1,2,3)) +
  scale_y_continuous(limits = c(0,1),
                     labels = scales::percent_format()) +
  theme_bw() +
  labs(x = NULL,
       y = "Estimated Percnet of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Party ID",
       shape = "Party ID")
## saved as figures/party_time.pdf
## top right panel of figure 2
## 5x7
## for caption: other variables held
## constant: ses 3, Boomer, White, male 
## points slightly offset to aid in inspection
## estimates referenced in paper
estimates.party
## 1986-1990 Republican: 68.5%, Democrats 87.1%
## 1991 - 1995 Republican 71.8%, 88.3%
## 1996-2000 Republican 60.8, 82.9%
## 2001-2005   Republican 51.6%, 87.0%
.718 - .883 ## 16.5 percentage point gap
.87 - .516 ## 35.4 percentage point gap

## non parametric estimate for comparison:
data.full %>%
  filter(year_5 == "1991-1995") %>%
  filter(party.harm %in% c("Republican",
                           "Democrat")) %>%
  group_by(party.harm) %>%
  summarize(mean = mean(inequality.num))
.876 - .683 ## 19.3 pp diff

data.full %>%
  filter(year_5 == "2001-2005") %>%
  filter(party.harm %in% c("Republican",
                           "Democrat")) %>%
  group_by(party.harm) %>%
  summarize(mean = mean(inequality.num))

.85 - .469 ## 38.1



## model 3: year*cohort

vcov <- vcov(full3a)
betas <- full3a$coefficients

set.seed(08540)
beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)

## create matrix
length(unique(data.cohortintx$year_5))
length(unique(data.cohortintx$cohort2))

matrix <- matrix(0, nrow = 24, 
                 ncol = length(beta.draws[1,]))
matrix[,1] <- 1
colnames(matrix) <- names(beta.draws[1,])
rownames(matrix) <- 1:nrow(matrix)

## year - different year for each level
## main year effects:

## for each level:
## each level gets a row of zeros (for first 1966-1970)
## and then a single year 
for(i in 2:6){
  matrix[i, i] <- 1
  matrix[i+6, i] <- 1
  matrix[i+12, i] <- 1
  matrix[i+18, i] <- 1
} 

## year*cohort interaction effects:
for(i in 2:6){
  matrix[i+6, i + 16] <- 1
  matrix[i+12, i + 21] <- 1
  matrix[i+18, i +26] <- 1
}

## main ses effects:
matrix[7:12, "cohort2Greatest"] <- 1
matrix[13:18, "cohort2Silent"] <- 1
matrix[19:24, "cohort2Gen-X"] <- 1

## other variables:
## set at Republican, ses3, White, male
matrix[,"ses3"] <- 1
matrix[,"race2Yes"] <- 1
matrix[, "gendermale"] <- 1

## for each draw of coefficients from the
## estimated distribution, we calculate
## estimated values 

sim <- base::apply(beta.draws, 1, function(x){
  vector <- matrix %*% x
  return(vector)
})
sim <- t(sim)

## calculate 2.5th and 97.5th percentile
## of estimated values for confidence intervals
estimates.cohort <- estimates.fct(sim)
estimates.cohort$frame <- "cohort"

data.cohortintx$year_5 <- droplevels(data.cohortintx$year_5)
data.cohortintx$cohort2 <- droplevels(data.cohortintx$cohort2)

estimates.cohort$year <- rep(levels(data.cohortintx$year_5), 4)
estimates.cohort$value <- c(rep(levels(data.cohortintx$cohort2)[1], 6), 
                            rep(levels(data.cohortintx$cohort2)[2], 6),
                            rep(levels(data.cohortintx$cohort2)[3], 6),
                            rep(levels(data.cohortintx$cohort2)[4], 6))

estimates.cohort$year_date <- gsub("-[0-9]{4}",
                                   "", 
                                   estimates.cohort$year)
estimates.cohort$year_date <- as.Date(paste0(estimates.cohort$year_date,
                                             "-01-01",
                                             sep = ""))

ggplot(estimates.cohort,
       aes(x = as.Date(year_date),
           y = Est,
           color = value,
           group = value,
           ymin = Lower,
           ymax = Upper)) +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  geom_line(position = position_dodge(2 * 365)) +
  scale_color_manual(values = colorblind_pal()(8)[c(1,2,3, 4)]) +
  scale_y_continuous(limits = c(0,1),
                     labels = scales::percent_format()) +
  scale_x_date(limits = c(as.Date("1966-01-01"),
                          as.Date("2011-12-30"))) +
  theme_bw() +
  labs(x = NULL,
       y = "Estimated Percnet of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Cohort")
## saved as figures/cohort_time.pdf
## lower left panel of figure 2
## 5x7
## for caption: other variables held
## constant: ses 3, Republican, White, male 
## points slightly offset to aid in inspection




## model 4: year*race

vcov <- vcov(full4a)
betas <- full4a$coefficients

set.seed(08540)
beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)

## create matrix

matrix <- matrix(0, nrow = 20, 
                 ncol = length(beta.draws[1,]))
matrix[,1] <- 1
colnames(matrix) <- names(beta.draws[1,])
rownames(matrix) <- 1:nrow(matrix)

## year - different year for each level
## main year effects:

## for each level:
## each level gets a row of zeros (for first 1966-1970)
## and then a single year 
for(i in 2:10){
  matrix[i, i] <- 1
  matrix[i+10, i] <- 1
} 

colnames(matrix)
## year*race interaction effects:
for(i in 2:10){
  matrix[i+10, i + 22] <- 1
}

## main race effects:
matrix[11:20, "race2Yes"] <- 1
matrix[,"ses3"] <- 1
matrix[, "gendermale"] <- 1
## other variables held at Republican,
## Boomers

## for each draw of coefficients from the
## estimated distribution, we calculate
## estimated values 

sim <- base::apply(beta.draws, 1, function(x){
  vector <- matrix %*% x
  return(vector)
})
sim <- t(sim)

## calculate 2.5th and 97.5th percentile
## of estimated values for confidence intervals
estimates.race <- estimates.fct(sim)
estimates.race$frame <- "race"

estimates.race$year <- rep(levels(data.full$year_5), 2)
estimates.race$value <- c(rep(levels(data.full$race2)[1], 10), 
                          rep(levels(data.full$race2)[2], 10))


estimates.race$year_date <- gsub("-[0-9]{4}",
                                 "", 
                                 estimates.race$year)
estimates.race$year_date <- as.Date(paste0(estimates.race$year_date,
                                           "-01-01",
                                           sep = ""))

estimates.race$value <- dplyr::recode(estimates.race$value,
                                      `No` = "Non-White",
                                      `Yes` = "White")

ggplot(estimates.race,
       aes(x = as.Date(year_date),
           y = Est,
           color = value,
           group = value,
           ymin = Lower,
           ymax = Upper)) +
  geom_point(position = position_dodge(2 * 365)) +
  geom_errorbar(width = 0,
                position = position_dodge(2 * 365)) +
  geom_line(position = position_dodge(2 * 365)) +
  scale_color_manual(values = colorblind_pal()(8)[c(6, 8)]) +
  scale_y_continuous(limits = c(0,1),
                     labels = scales::percent_format()) +
  theme_bw() +
  labs(x = NULL,
       y = "Estimated Percnet of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Race")
## saved as figures/race_time.pdf
## lower right panel of figure 2
## 5x7
## for caption: other variables held
## constant: ses 3, Republican, Boomers, male
## points slightly offset to aid in inspection


######################################
## ses-party-year interaction ########
######################################

## We create Figure 3 for main text
## in this block of code. Figure 3 visualizes
## the estimated percent of Americans who agree
## the "rich are getting richer" by party, SES,
## and 5-year bin 

full5 <- lm(inequality.num ~ year_5*party.harm*ses + race2 + cohort2 +
              gender, data = data.full)

vcov <- vcov(full5)
betas <- full5$coefficients

set.seed(08540)
beta.draws <- MASS::mvrnorm(10000, mu = betas, Sigma = vcov)

## create matrix 
matrix <- matrix(0, nrow = 150, 
                 ncol = length(beta.draws[1,]))
matrix[,1] <- 1
colnames(matrix) <- names(beta.draws[1,])
rownames(matrix) <- 1:nrow(matrix)

## year - different year for each level
## main year effects:

## for each level:
## each level gets a row of zeros (for first 1966-1970)
## and then a single year 
for(i in 2:10){
  matrix[i, i] <- 1
  matrix[i+10, i] <- 1
  matrix[i+20, i] <- 1
  matrix[i+30, i] <- 1
  matrix[i+40, i] <- 1
  matrix[i+50, i] <- 1
  matrix[i+60, i] <- 1
  matrix[i+70, i] <- 1
  matrix[i+80, i] <- 1
  matrix[i+90, i] <- 1
  matrix[i+100, i] <- 1
  matrix[i+110, i] <- 1
  matrix[i+120, i] <- 1
  matrix[i+130, i] <- 1
  matrix[i+140, i] <- 1
} 

## party 
## first 50: republicans - blank
## second 50: democrats
## third 50: other
matrix[51:100, "party.harmDemocrat"] <- 1
matrix[101:150, "party.harmIndep/Other"] <- 1

## ses
## ## 1-10: ses 1 - blank
## 10-20: ses 2
## 20-30: ses 3
## 30-40: ses 4
## 40-50: ses 5
## repeat for each party group 
matrix[c(11:20,
         61:70,
         111:120), "ses2"] <- 1
matrix[c(21:30,
         71:80,
         121:130), "ses3"] <- 1
matrix[c(31:40,
         81:90,
         131:140), "ses4"] <- 1
matrix[c(41:50,
         91:100,
         141:150), "ses5"] <- 1


## ses:party interaction
## first 50 blank - republicans
matrix[61:70, "party.harmDemocrat:ses2"] <- 1
matrix[71:80, "party.harmDemocrat:ses3"] <- 1
matrix[81:90, "party.harmDemocrat:ses4"] <- 1
matrix[91:100, "party.harmDemocrat:ses5"] <- 1
matrix[111:120, "party.harmIndep/Other:ses2"] <- 1
matrix[121:130, "party.harmIndep/Other:ses3"] <- 1
matrix[131:140, "party.harmIndep/Other:ses4"] <- 1
matrix[141:150, "party.harmIndep/Other:ses5"] <- 1

## other variables set constant
matrix[,"race2Yes"] <- 1
matrix[, "gendermale"] <- 1
## Boomer (cohort2) 

## year*party interaction effects:
## first 50: republicans - blank
## second 50: democrats
## third 50: other
for(i in 2:10){
  matrix[i+50, i + 22] <- 1
  matrix[i+60, i + 22] <- 1
  matrix[i+70, i + 22] <- 1
  matrix[i+80, i + 22] <- 1
  matrix[i+90, i + 22] <- 1
  matrix[i+100, i + 31] <- 1
  matrix[i+110, i + 31] <- 1
  matrix[i+120, i + 31] <- 1
  matrix[i+130, i + 31] <- 1
  matrix[i+140, i + 31] <- 1
}


## year*ses interaction effects:
## 1-10: ses 1 - blank
## 10-20: ses 2
## 20-30: ses 3
## 30-40: ses 4
## 40-50: ses 5
## repeat for each party group 

for(i in 2:10){
  matrix[i+10, i + 40] <- 1
  matrix[i+20, i + 49] <- 1
  matrix[i+30, i + 58] <- 1
  matrix[i+40, i + 67] <- 1
  matrix[i+60, i + 40] <- 1
  matrix[i+70, i + 49] <- 1
  matrix[i+80, i + 58] <- 1
  matrix[i+90, i + 67] <- 1
  matrix[i+110, i + 40] <- 1
  matrix[i+120, i + 49] <- 1
  matrix[i+130, i + 58] <- 1
  matrix[i+140, i + 67] <- 1
}


## party-year-ses effects
matrix[61,"year_51971-1975:party.harmDemocrat:ses2"] <- 1
matrix[62,"year_51976-1980:party.harmDemocrat:ses2"] <- 1
matrix[63,"year_51981-1985:party.harmDemocrat:ses2"] <- 1
matrix[64,"year_51986-1990:party.harmDemocrat:ses2"] <- 1
matrix[65,"year_51991-1995:party.harmDemocrat:ses2"] <- 1
matrix[66,"year_51996-2000:party.harmDemocrat:ses2"] <- 1
matrix[67,"year_51996-2000:party.harmDemocrat:ses2"] <- 1
matrix[68,"year_52001-2005:party.harmDemocrat:ses2"] <- 1
matrix[69,"year_52006-2010:party.harmDemocrat:ses2"] <- 1
matrix[70,"year_52011-2013:party.harmDemocrat:ses2"] <- 1
matrix[71,"year_51971-1975:party.harmDemocrat:ses3"] <- 1
matrix[72,"year_51976-1980:party.harmDemocrat:ses3"] <- 1
matrix[73,"year_51981-1985:party.harmDemocrat:ses3"] <- 1
matrix[74,"year_51986-1990:party.harmDemocrat:ses3"] <- 1
matrix[75,"year_51991-1995:party.harmDemocrat:ses3"] <- 1
matrix[76,"year_51996-2000:party.harmDemocrat:ses3"] <- 1
matrix[77,"year_51996-2000:party.harmDemocrat:ses3"] <- 1
matrix[78,"year_52001-2005:party.harmDemocrat:ses3"] <- 1
matrix[79,"year_52006-2010:party.harmDemocrat:ses3"] <- 1
matrix[80,"year_52011-2013:party.harmDemocrat:ses3"] <- 1
matrix[81,"year_51971-1975:party.harmDemocrat:ses4"] <- 1
matrix[82,"year_51976-1980:party.harmDemocrat:ses4"] <- 1
matrix[83,"year_51981-1985:party.harmDemocrat:ses4"] <- 1
matrix[84,"year_51986-1990:party.harmDemocrat:ses4"] <- 1
matrix[85,"year_51991-1995:party.harmDemocrat:ses4"] <- 1
matrix[86,"year_51996-2000:party.harmDemocrat:ses4"] <- 1
matrix[87,"year_51996-2000:party.harmDemocrat:ses4"] <- 1
matrix[88,"year_52001-2005:party.harmDemocrat:ses4"] <- 1
matrix[89,"year_52006-2010:party.harmDemocrat:ses4"] <- 1
matrix[90,"year_52011-2013:party.harmDemocrat:ses4"] <- 1
matrix[91,"year_51971-1975:party.harmDemocrat:ses5"] <- 1
matrix[92,"year_51976-1980:party.harmDemocrat:ses5"] <- 1
matrix[93,"year_51981-1985:party.harmDemocrat:ses5"] <- 1
matrix[94,"year_51986-1990:party.harmDemocrat:ses5"] <- 1
matrix[95,"year_51991-1995:party.harmDemocrat:ses5"] <- 1
matrix[96,"year_51996-2000:party.harmDemocrat:ses5"] <- 1
matrix[97,"year_51996-2000:party.harmDemocrat:ses5"] <- 1
matrix[98,"year_52001-2005:party.harmDemocrat:ses5"] <- 1
matrix[99,"year_52006-2010:party.harmDemocrat:ses5"] <- 1
matrix[100,"year_52011-2013:party.harmDemocrat:ses5"] <- 1
matrix[111,"year_51971-1975:party.harmIndep/Other:ses2"] <- 1
matrix[112,"year_51976-1980:party.harmIndep/Other:ses2"] <- 1
matrix[113,"year_51981-1985:party.harmIndep/Other:ses2"] <- 1
matrix[114,"year_51986-1990:party.harmIndep/Other:ses2"] <- 1
matrix[115,"year_51991-1995:party.harmIndep/Other:ses2"] <- 1
matrix[116,"year_51996-2000:party.harmIndep/Other:ses2"] <- 1
matrix[117,"year_51996-2000:party.harmIndep/Other:ses2"] <- 1
matrix[118,"year_52001-2005:party.harmIndep/Other:ses2"] <- 1
matrix[119,"year_52006-2010:party.harmIndep/Other:ses2"] <- 1
matrix[120,"year_52011-2013:party.harmIndep/Other:ses2"] <- 1
matrix[121,"year_51971-1975:party.harmIndep/Other:ses3"] <- 1
matrix[122,"year_51976-1980:party.harmIndep/Other:ses3"] <- 1
matrix[123,"year_51981-1985:party.harmIndep/Other:ses3"] <- 1
matrix[124,"year_51986-1990:party.harmIndep/Other:ses3"] <- 1
matrix[125,"year_51991-1995:party.harmIndep/Other:ses3"] <- 1
matrix[126,"year_51996-2000:party.harmIndep/Other:ses3"] <- 1
matrix[127,"year_51996-2000:party.harmIndep/Other:ses3"] <- 1
matrix[128,"year_52001-2005:party.harmIndep/Other:ses3"] <- 1
matrix[129,"year_52006-2010:party.harmIndep/Other:ses3"] <- 1
matrix[130,"year_52011-2013:party.harmIndep/Other:ses3"] <- 1
matrix[131,"year_51971-1975:party.harmIndep/Other:ses4"] <- 1
matrix[132,"year_51976-1980:party.harmIndep/Other:ses4"] <- 1
matrix[133,"year_51981-1985:party.harmIndep/Other:ses4"] <- 1
matrix[134,"year_51986-1990:party.harmIndep/Other:ses4"] <- 1
matrix[135,"year_51991-1995:party.harmIndep/Other:ses4"] <- 1
matrix[136,"year_51996-2000:party.harmIndep/Other:ses4"] <- 1
matrix[137,"year_51996-2000:party.harmIndep/Other:ses4"] <- 1
matrix[138,"year_52001-2005:party.harmIndep/Other:ses4"] <- 1
matrix[139,"year_52006-2010:party.harmIndep/Other:ses4"] <- 1
matrix[140,"year_52011-2013:party.harmIndep/Other:ses4"] <- 1
matrix[141,"year_51971-1975:party.harmIndep/Other:ses5"] <- 1
matrix[142,"year_51976-1980:party.harmIndep/Other:ses5"] <- 1
matrix[143,"year_51981-1985:party.harmIndep/Other:ses5"] <- 1
matrix[144,"year_51986-1990:party.harmIndep/Other:ses5"] <- 1
matrix[145,"year_51991-1995:party.harmIndep/Other:ses5"] <- 1
matrix[146,"year_51996-2000:party.harmIndep/Other:ses5"] <- 1
matrix[147,"year_51996-2000:party.harmIndep/Other:ses5"] <- 1
matrix[148,"year_52001-2005:party.harmIndep/Other:ses5"] <- 1
matrix[149,"year_52006-2010:party.harmIndep/Other:ses5"] <- 1
matrix[150,"year_52011-2013:party.harmIndep/Other:ses5"] <- 1


## for each draw of coefficients from the
## estimated distribution, we calculate
## estimated values 

sim <- base::apply(beta.draws, 1, function(x){
  vector <- matrix %*% x
  return(vector)
})
sim <- t(sim)

## calculate 2.5th and 97.5th percentile
## of estimated values for confidence intervals
estimates.intx <- estimates.fct(sim)

estimates.intx$year <- rep(levels(data.full$year_5), 15)
estimates.intx$party <- c(rep(levels(data.full$party.harm)[1], 50), 
                          rep(levels(data.full$party.harm)[2], 50),
                          rep(levels(data.full$party.harm)[3], 50))
estimates.intx$ses <- c(rep(levels(data.full$ses)[1], 10), 
                        rep(levels(data.full$ses)[2], 10),
                        rep(levels(data.full$ses)[3], 10),
                        rep(levels(data.full$ses)[4], 10),
                        rep(levels(data.full$ses)[5], 10),
                        rep(levels(data.full$ses)[1], 10), 
                        rep(levels(data.full$ses)[2], 10),
                        rep(levels(data.full$ses)[3], 10),
                        rep(levels(data.full$ses)[4], 10),
                        rep(levels(data.full$ses)[5], 10),
                        rep(levels(data.full$ses)[1], 10), 
                        rep(levels(data.full$ses)[2], 10),
                        rep(levels(data.full$ses)[3], 10),
                        rep(levels(data.full$ses)[4], 10),
                        rep(levels(data.full$ses)[5], 10))



estimates.intx$year_date <- gsub("-[0-9]{4}",
                                 "", 
                                 estimates.intx$year)
estimates.intx$year_date <- as.Date(paste0(estimates.intx$year_date,
                                           "-01-01",
                                           sep = ""))



ggplot(estimates.intx,
       aes(x = year_date,
           y = Est,
           ymin = Lower,
           ymax = Upper,
           shape = party,
           color = party)) +
  geom_point(position = position_dodge(width = 2*365)) +
  geom_line(position = position_dodge(width = 2*365)) +
  geom_errorbar(width = 0,
                position = position_dodge(width = 2*365)) +
  facet_wrap(~ ses, ncol = 2) +
  theme_bw() +
  labs(x = "Minimum Year in 5 Year Range",
       y = "Estimated Percent of Respondents\nWho Think 'Rich Getting Richer'",
       color = "Party",
       shape = "Party") +
  scale_color_manual(values = c("#2c7fb8",
                                "#31a354",
                                "#f03b20")) +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = c(0,.25,.5,.75,1),
                     limits = c(0,1))
## saved as figures/ses_party_intx
## Figure 3 in main text
## 8x7

## Black and White version: 
ggplot(estimates.intx,
       aes(x = year_date,
           y = Est,
           ymin = Lower,
           ymax = Upper,
           shape = party,
           lty = party)) +
  geom_point(position = position_dodge(width = 2*365)) +
  geom_line(position = position_dodge(width = 2*365)) +
  geom_errorbar(width = 0,
                position = position_dodge(width = 2*365)) +
  facet_wrap(~ ses, ncol = 2) +
  theme_bw() +
  labs(x = "Minimum Year in 5 Year Range",
       y = "Estimated Percent of Respondents\nWho Think 'Rich Getting Richer'",
       shape = "Party",
       lty = "Party") +
  theme(legend.position = "bottom") +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = c(0,.25,.5,.75,1),
                     limits = c(0,1))
## figures/ses_party_intx_bw
## Figure 3 in main text
## 8x7
